Meta-analysis of the Interoceptive Accuracy Scale

Study 1

library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Data Preparation

Sample 1

Murphy (https://osf.io/3m5nh)

df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

Sample 2

Gaggero (2021) (https://osf.io/5x9sg)

load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)

Sample 3

Campos (2022) - Study 1 (https://osf.io/j6ef3)

df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 4

Todd (2022) - https://osf.io/ms354/

df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

Sample 5

Arslanova (2022) - https://osf.io/mp3cy/

df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |> 
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  select(-starts_with("Participant")) |> 
  mutate_all(as.numeric)

Sample 6

Brand (2022) - https://osf.io/xwz6g/

df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 7

Brand (2023) - https://osf.io/3f2h6/

df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
  )

Sample 8

Lin (2023) - https://osf.io/3eztd - Note: No tickle because same Chinese character

df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH)

Sample 9

VonMohr (2023) - Study 3 (https://osf.io/7p9u5/)

df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

Sample 10

Makowski

• Note: First sample has some missing items
• Note: Analog scales

df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b)
save(data, file = "../data/data.RData")

Participants

• Sample 1a: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
• Sample 1b: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
• TODO: ADD INFO

Total N = 31317.

Descriptive

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10a",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 10b",
         x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))


p1 <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b) |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b")) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  scale_color_material_d() +
  scale_linetype_manual(values = c("solid", "dashed"), guide="none") +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        plot.title = element_text(face="bold")) +
  labs(title = "Distribution of responses for all items across various datasets", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free_y", nrow = 5)

ggsave("figures/Figure1.png", p1, width=10, height=10, dpi=200, bg="white")
p1

data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))

data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), mutate(data8b, Tickle = NA), data9,
  mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA), data10b
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

All data

make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

make_correlation(data_all)

Sample 1a

make_correlation(data1a)

Sample 1b

make_correlation(data1b)

Sample 2

make_correlation(data2)

Sample 3

make_correlation(data3)

Sample 4

make_correlation(data4)

Sample 5

make_correlation(data5)

Sample 6

make_correlation(data6)

Sample 7a

make_correlation(data7a)

Sample 7b

make_correlation(data7b)

Sample 7c

make_correlation(data7c)

Sample 8a

make_correlation(data8a)

Sample 8b

make_correlation(data8b)

Sample 9

make_correlation(data9)

Sample 10a

make_correlation(data10a)

Sample 10b

make_correlation(data10b)

Graph Analysis

Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

All data

uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.368

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.293
 Urinate Defecate 0.270

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.234
  Heart Breathing 0.225
 Hungry   Thirsty 0.218

uva$keep_remove

$keep
[1] "Itch"

$remove
[1] "Tickle"

Sample 1a

uva <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241

uva$keep_remove

NULL

Sample 1b

uva <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"

Sample 3

uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208

uva$keep_remove

$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"

Sample 4

uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.415
 Tickle      Itch 0.351

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.293
 Urinate Defecate 0.282
    Wind     Burp 0.275
  Sneeze    Cough 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i node_j   wto
 Muscles   Pain 0.205

uva$keep_remove

$keep
[1] "Heart" "Itch" 

$remove
[1] "Breathing" "Tickle"   

Sample 5

uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.251

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
 Bruise   Itch 0.241

uva$keep_remove

NULL

Sample 6

uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237

uva$keep_remove

$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  

Sample 7a

uva <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 

Sample 7b

uva <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237

uva$keep_remove

$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 

Sample 7c

uva <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218

uva$keep_remove

$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"

Sample 8a

uva <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213

uva$keep_remove

NULL

Sample 8b

uva <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206

uva$keep_remove

NULL

Sample 9

uva <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24

uva$keep_remove

$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 

Sample 10a

uva <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209

uva$keep_remove

NULL

Sample 10b

uva <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

uva$keep_remove

NULL

Networks

All data

make_egas <- function(data) {
  egas <- list()
  for (model in c("glasso")) {  # , "TMFG"
    for (algo in c("walktrap", "louvain")) {
      for (type in c("ega", "ega.fit")) { # "" # "hierega"
        if (type == "ega.fit" & algo == "louvain") next # Too slow
        egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
          data = data,
          seed = 123,
          model = model,
          algorithm = algo,
          EGA.type = type,
          type = "resampling",
          plot.itemStability = FALSE,
          verbose = FALSE
        )
      }
    }
  }

  p <- EGAnet::compare.EGA.plots(
    egas$glasso_walktrap_ega, 
    # egas$glasso_walktrap_ega.fit,
    egas$glasso_louvain_ega, 
    # egas$TMFG_louvain_ega,
    # egas$TMFG_walktrap_ega, 
    # egas$TMFG_walktrap_ega.fit,
    labels = c(
      "glasso_walktrap_ega", 
      # "glasso_walktrap_ega.fit",
      "glasso_louvain_ega" 
      # "TMFG_louvain_ega",
      # "TMFG_walktrap_ega", 
      # "TMFG_walktrap_ega.fit"
    ),
    rows = 3,
    plot.all = FALSE
  )$all
  list(egas = egas, p = p)
}


egas_all <- make_egas(data_all)
egas_all$p

Sample 1a

egas1a <- make_egas(data1a)
egas1a$p

Sample 1b

egas1b <- make_egas(data1b)
egas1b$p

Sample 2

egas2 <- make_egas(data2)
egas2$p

Sample 3

egas3 <- make_egas(data3)
egas3$p

Sample 4

egas4 <- make_egas(data4)
egas4$p

Sample 5

egas5 <- make_egas(data5)
egas5$p

Sample 6

egas6 <- make_egas(data6)
egas6$p

Sample 7a

egas7a <- make_egas(data7a)
egas7a$p

Sample 7b

egas7b <- make_egas(data7b)
egas7b$p

Sample 7c

egas7c <- make_egas(data7c)
egas7c$p

Sample 8a

egas8a <- make_egas(data8a)
egas8a$p

Sample 8b

egas8b <- make_egas(data8b)
egas8b$p

Sample 9

egas9 <- make_egas(data9)
egas9$p

Sample 10a

egas10a <- make_egas(data10a)
egas10a$p

Sample 10b

egas10b <- make_egas(data10b)
egas10b$p

Structure Stability

All data

patchwork::wrap_plots(lapply(egas_all$egas, plot), nrow = 3)

Sample 1a

patchwork::wrap_plots(lapply(egas1a$egas, plot), nrow = 3)

Sample 1b

patchwork::wrap_plots(lapply(egas1b$egas, plot), nrow = 3)

Sample 2

patchwork::wrap_plots(lapply(egas2$egas, plot), nrow = 3)

Sample 3

patchwork::wrap_plots(lapply(egas3$egas, plot), nrow = 3)

Sample 4

patchwork::wrap_plots(lapply(egas4$egas, plot), nrow = 3)

Sample 5

patchwork::wrap_plots(lapply(egas5$egas, plot), nrow = 3)

Sample 6

patchwork::wrap_plots(lapply(egas6$egas, plot), nrow = 3)

Sample 7a

patchwork::wrap_plots(lapply(egas7a$egas, plot), nrow = 3)

Sample 7b

patchwork::wrap_plots(lapply(egas7b$egas, plot), nrow = 3)

Sample 7c

patchwork::wrap_plots(lapply(egas7c$egas, plot), nrow = 3)

Sample 8a

patchwork::wrap_plots(lapply(egas8a$egas, plot), nrow = 3)

Sample 8b

patchwork::wrap_plots(lapply(egas8b$egas, plot), nrow = 3)

Sample 9

patchwork::wrap_plots(lapply(egas9$egas, plot), nrow = 3)

Sample 10a

patchwork::wrap_plots(lapply(egas10a$egas, plot), nrow = 3)

Sample 10b

patchwork::wrap_plots(lapply(egas10b$egas, plot), nrow = 3)

Scores

ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)

get_scores <- function(data, ega) {
  EGAnet::net.scores(data, ega)$scores$std.scores |>
    as.data.frame() |>
    setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))
}

ega_scores_1a <- get_scores(data1a, ega)
ega_scores_1b <- get_scores(data1b, ega)
ega_scores_2 <- get_scores(data2, ega)
ega_scores_3 <- get_scores(data3, ega)
ega_scores_4 <- get_scores(data4, ega)
ega_scores_5 <- get_scores(data5, ega)
ega_scores_6 <- get_scores(data6, ega)
ega_scores_7a <- get_scores(data7a, ega)
ega_scores_7b <- get_scores(data7b, ega)
ega_scores_7c <- get_scores(data7c, ega)
ega_scores_8a <- get_scores(data8a, ega)
ega_scores_8b <- get_scores(data8b, ega)
ega_scores_9 <- get_scores(data9, ega)
ega_scores_10a <- get_scores(data10a, ega)
ega_scores_10b <- get_scores(data10b, ega)

Factor Analysis

How Many Factors

All data

n <- parameters::n_factors(data_all, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 1a

n <- parameters::n_factors(data1a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 1b

n <- parameters::n_factors(data1b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
5	BIC (adjusted)	BIC

Sample 2

n <- parameters::n_factors(data2, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
4	BIC	BIC
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree

Sample 3

n <- parameters::n_factors(data3, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
3	CNG	CNG
3	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
4	BIC	BIC

Sample 4

n <- parameters::n_factors(data4, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
5	BIC	BIC

Sample 5

n <- parameters::n_factors(data5, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Optimal coordinates	Scree
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
1	BIC	BIC
3	CNG	CNG
3	BIC (adjusted)	BIC
4	Bartlett	Barlett
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Anderson	Barlett
5	Lawley	Barlett

Sample 6

n <- parameters::n_factors(data6, n_max = 12)

Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected.  Examine the results carefully

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 7a

n <- parameters::n_factors(data7a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
3	BIC	BIC
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 7b

n <- parameters::n_factors(data7b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Bentler	Bentler
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree

Sample 7c

n <- parameters::n_factors(data7c, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Scree (SE)	Scree_SE

Sample 8a

n <- parameters::n_factors(data8a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	Bentler	Bentler
3	CNG	CNG
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree
5	BIC	BIC

Sample 8b

n <- parameters::n_factors(data8b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	BIC	BIC
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree
5	Scree (SE)	Scree_SE
5	BIC (adjusted)	BIC

Sample 9

n <- parameters::n_factors(data9, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 10a

n <- parameters::n_factors(data10a, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
3	Scree (SE)	Scree_SE
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	BIC	BIC
4	BIC (adjusted)	BIC

Sample 10b

n <- parameters::n_factors(data10b, n_max = 12)

plot(n)

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
3	CNG	CNG
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
5	BIC (adjusted)	BIC

EFA

All data

efa5 <- parameters::factor_analysis(data_all, n = 4, rotation = "oblimin", sort = TRUE)

Loading required namespace: GPArotation

plot(efa5)

display(efa5)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Burp	0.74	0.03	-0.02	-0.07	1.02	0.48
Cough	0.70	-0.02	0.02	0.06	1.02	0.47
Wind	0.67	0.02	-0.02	-2.95e-03	1.00	0.55
Sneeze	0.63	-0.06	0.04	0.14	1.12	0.52
Vomit	0.44	0.06	0.04	0.18	1.39	0.63
Temp	0.27	0.24	0.05	0.22	2.98	0.61
Sex_arousal	0.26	0.22	0.03	0.20	2.90	0.67
Taste	0.24	0.20	0.10	0.15	2.99	0.70
Breathing	0.05	0.61	-0.04	0.04	1.03	0.59
Hungry	-0.11	0.54	-5.62e-03	0.19	1.34	0.64
Heart	0.08	0.49	0.03	-0.07	1.10	0.73
Thirsty	-0.07	0.47	-7.04e-03	0.25	1.60	0.64
Pain	0.16	0.40	0.10	0.10	1.60	0.62
Muscles	0.30	0.36	0.07	0.03	2.03	0.59
Blood_Sugar	0.11	0.35	0.22	-0.18	2.47	0.76
Bruise	0.25	0.34	0.22	-0.19	3.30	0.65
Tickle	-0.04	-0.04	0.82	0.04	1.02	0.37
Itch	0.02	6.24e-03	0.75	-0.04	1.01	0.43
Affective_touch	0.08	0.19	0.33	0.11	2.00	0.70
Urinate	0.04	0.10	0.04	0.65	1.06	0.45
Defecate	0.15	0.04	0.05	0.61	1.14	0.47

The 4 latent factors (oblimin rotation) accounted for 41.65% of the total variance of the original data (MR1 = 14.63%, MR3 = 11.59%, MR2 = 8.13%, MR4 = 7.30%).

Sample 1a

efa4 <- parameters::factor_analysis(data1a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR2	MR3	MR1	MR4	Complexity	Uniqueness
Tickle	0.72	0.03	-0.10	0.11	1.09	0.47
Itch	0.67	0.14	-9.44e-03	-0.10	1.13	0.46
Affective_touch	0.49	-0.13	0.29	0.13	1.97	0.59
Blood_Sugar	0.29	0.21	0.18	-0.23	3.51	0.75
Burp	-0.03	0.73	-0.03	-1.21e-03	1.00	0.50
Cough	7.56e-03	0.69	7.00e-03	-0.01	1.00	0.52
Sneeze	0.07	0.54	7.09e-03	0.20	1.30	0.55
Wind	0.08	0.46	0.07	0.16	1.37	0.62
Hungry	-0.16	0.04	0.65	0.03	1.13	0.60
Thirsty	-4.84e-03	0.01	0.49	0.11	1.11	0.70
Breathing	0.14	0.05	0.44	0.03	1.24	0.70
Heart	0.11	-0.07	0.40	0.08	1.32	0.78
Bruise	0.31	0.25	0.34	-0.25	3.74	0.59
Pain	0.29	0.01	0.32	0.20	2.67	0.61
Temp	0.19	0.13	0.29	0.16	2.91	0.67
Muscles	0.20	0.13	0.26	0.16	3.24	0.69
Sex_arousal	0.18	0.02	0.21	0.15	2.81	0.82
Defecate	0.10	0.11	2.77e-03	0.66	1.10	0.46
Urinate	-0.11	0.13	0.24	0.52	1.65	0.54
Taste	0.23	0.06	0.17	0.30	2.60	0.69
Vomit	0.21	0.18	0.11	0.25	3.26	0.71

The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).

Sample 1b

efa4 <- parameters::factor_analysis(data1b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR3	MR4	Complexity	Uniqueness
Urinate	0.69	-0.18	0.13	-0.04	1.22	0.54
Defecate	0.68	-0.13	0.04	0.06	1.10	0.55
Breathing	0.56	0.03	0.06	7.96e-03	1.03	0.63
Hungry	0.55	0.15	-7.15e-03	0.02	1.16	0.60
Thirsty	0.52	0.05	0.14	-0.01	1.16	0.62
Sex_arousal	0.50	0.14	-0.04	0.20	1.49	0.56
Temp	0.46	0.23	0.16	-0.01	1.77	0.52
Pain	0.46	0.43	-0.05	-0.03	2.03	0.49
Taste	0.40	0.21	-0.07	0.12	1.77	0.69
Heart	0.40	0.06	-0.04	0.15	1.36	0.76
Muscles	0.35	0.24	0.07	0.09	2.07	0.65
Vomit	0.32	0.06	0.23	0.11	2.20	0.68
Tickle	0.02	0.69	0.01	0.04	1.01	0.47
Itch	-0.02	0.67	0.18	0.03	1.16	0.41
Bruise	0.07	0.44	0.15	0.06	1.32	0.65
Blood_Sugar	-0.08	0.38	0.16	0.09	1.61	0.77
Affective_touch	0.38	0.38	-0.13	0.05	2.26	0.63
Sneeze	0.05	0.05	0.71	0.01	1.02	0.41
Cough	0.04	0.02	0.71	0.04	1.01	0.43
Wind	7.52e-03	-0.04	-5.18e-03	0.94	1.00	0.15
Burp	-0.01	0.27	0.17	0.47	1.91	0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

CFA

We discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant.

Ambiguous: - Temp - Vomit - Affective_touch - Sex_arousal - Taste

Lower-Level Factors

All data

m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"

m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"

m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"

is_converged <- function(m) {
  ev <- eigen(lavaan::lavTech(m, "cov.lv")[[1]], symmetric=TRUE, only.values=TRUE)$values
  if (any(ev < -1 * .Machine$double.eps^(3/4))) {
    "Not converged"
  } else {
    ""
  }
}

compare_mods <- function(data) {
  m1 <- lavaan::cfa(m1, data = data)
  m4 <- lavaan::cfa(m4, data = data)
  m5 <- lavaan::cfa(m5, data = data)
  m6 <- lavaan::cfa(m6, data = data)
  m7 <- lavaan::cfa(m7, data = data)
  
  anova(m1, m4, m5, m6, m7) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m1, m4, m5, m6, m7), is_converged))
}

display(compare_mods(data_all))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-168983	-168574	2080.93	56
m6	-164335	-163977	6740.02	62	< .001
m5	-160047	-159731	11038.16	67	< .001
m4	-154973	-154690	16120.58	71	< .001
m1	-141854	-141621	29251.50	77	< .001

m7_all <- lavaan::cfa(m7, data = data_all)
parameters::parameters(m7_all, standardize=TRUE) |>
  display()

# Loading

Link	Coefficient	SE	95% CI	z	p
HungryThirsty =~ Hungry	0.68	4.79e-03	(0.67, 0.69)	142.07	< .001
HungryThirsty =~ Thirsty	0.71	4.77e-03	(0.70, 0.72)	148.26	< .001
UrinateDefecate =~ Urinate	0.77	3.94e-03	(0.76, 0.78)	195.43	< .001
UrinateDefecate =~ Defecate	0.77	3.93e-03	(0.77, 0.78)	196.53	< .001
ItchBruise =~ Itch	0.50	5.73e-03	(0.49, 0.51)	87.60	< .001
ItchBruise =~ Bruise	0.67	5.99e-03	(0.66, 0.68)	112.01	< .001
MusclesPain =~ Muscles	0.71	4.15e-03	(0.70, 0.72)	171.30	< .001
MusclesPain =~ Pain	0.68	4.22e-03	(0.68, 0.69)	161.84	< .001
HeartBreathing =~ Breathing	0.79	5.19e-03	(0.78, 0.80)	151.51	< .001
HeartBreathing =~ Heart	0.59	5.13e-03	(0.58, 0.60)	114.83	< .001
CoughSneeze =~ Cough	0.81	3.49e-03	(0.80, 0.81)	230.73	< .001
CoughSneeze =~ Sneeze	0.75	3.65e-03	(0.74, 0.76)	205.71	< .001
WindBurp =~ Wind	0.77	3.58e-03	(0.76, 0.78)	215.22	< .001
WindBurp =~ Burp	0.82	3.44e-03	(0.81, 0.83)	238.24	< .001

# Correlation

Link	Coefficient	SE	95% CI	z	p
HungryThirsty ~~ UrinateDefecate	0.65	6.26e-03	(0.63, 0.66)	103.09	< .001
HungryThirsty ~~ ItchBruise	0.45	9.03e-03	(0.43, 0.46)	49.41	< .001
HungryThirsty ~~ MusclesPain	0.63	7.04e-03	(0.61, 0.64)	89.04	< .001
HungryThirsty ~~ HeartBreathing	0.64	7.01e-03	(0.63, 0.66)	91.55	< .001
HungryThirsty ~~ CoughSneeze	0.45	7.14e-03	(0.44, 0.47)	63.70	< .001
HungryThirsty ~~ WindBurp	0.43	7.12e-03	(0.42, 0.44)	60.47	< .001
UrinateDefecate ~~ ItchBruise	0.43	8.38e-03	(0.42, 0.45)	51.81	< .001
UrinateDefecate ~~ MusclesPain	0.64	6.23e-03	(0.63, 0.66)	103.29	< .001
UrinateDefecate ~~ HeartBreathing	0.53	6.81e-03	(0.52, 0.55)	78.24	< .001
UrinateDefecate ~~ CoughSneeze	0.56	6.03e-03	(0.54, 0.57)	92.16	< .001
UrinateDefecate ~~ WindBurp	0.50	6.19e-03	(0.49, 0.52)	81.27	< .001
ItchBruise ~~ MusclesPain	0.80	7.96e-03	(0.79, 0.82)	100.88	< .001
ItchBruise ~~ HeartBreathing	0.51	8.74e-03	(0.49, 0.53)	58.21	< .001
ItchBruise ~~ CoughSneeze	0.64	7.69e-03	(0.62, 0.65)	82.74	< .001
ItchBruise ~~ WindBurp	0.64	7.55e-03	(0.62, 0.65)	84.63	< .001
MusclesPain ~~ HeartBreathing	0.64	6.99e-03	(0.62, 0.65)	91.42	< .001
MusclesPain ~~ CoughSneeze	0.66	6.06e-03	(0.65, 0.68)	109.42	< .001
MusclesPain ~~ WindBurp	0.63	6.10e-03	(0.62, 0.65)	103.72	< .001
HeartBreathing ~~ CoughSneeze	0.52	6.81e-03	(0.50, 0.53)	75.70	< .001
HeartBreathing ~~ WindBurp	0.45	6.95e-03	(0.44, 0.47)	65.21	< .001
CoughSneeze ~~ WindBurp	0.72	4.84e-03	(0.71, 0.73)	149.24	< .001

Sample 1a

display(compare_mods(data1a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1951	-1749	88.22	56
m6	-1925	-1748	126.29	62	< .001
m5	-1860	-1703	201.38	67	< .001	Not converged
m4	-1834	-1694	234.82	71	< .001
m1	-1661	-1546	420.24	77	< .001

Conclusion: The 7th factor solution was preferred in most datasets, including by indices that penalize increased number of parameters (such as BIC).

Higher Order Factors

All data

m7h1 <- paste0(m7, "
              Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
              ")

m7h2 <- paste0(m7, "
              Expulsion =~ CoughSneeze + WindBurp
              Interoception =~ HungryThirsty + UrinateDefecate + HeartBreathing
              ")

m7h3 <- paste0(m7, "
              Valenced =~ MusclesPain + HeartBreathing + ItchBruise
              Expulsion =~ CoughSneeze + WindBurp
              Homeostasis =~ HungryThirsty + UrinateDefecate
              ")

compare_mods_h <- function(data) {
  m7 <- lavaan::cfa(m7, data = data)
  m7h2 <- lavaan::cfa(m7h2, data = data)
  m7h1 <- lavaan::cfa(m7h1, data = data)
  m7h3 <- lavaan::cfa(m7h3, data = data)
  
  anova(m7, m7h1, m7h2, m7h3) |>
    parameters::parameters() |>
    mutate(AIC = format(round(AIC, 0), nsmall = 0),
           BIC = format(round(BIC, 0), nsmall = 0),
           Chi2 = format(round(Chi2, 2), nsmall = 2),
           Converged = sapply(list(m7, m7h1, m7h2, m7h3), is_converged))
}

display(compare_mods_h(data_all))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-168983	-168574	2080.93	56
m7h2	-168224	-167899	2859.83	66	< .001
m7h3	-166932	-166615	4153.51	67	< .001
m7h1	-164501	-164210	6590.01	70	< .001

No evidence for higher order factors.

Sample 1a

display(compare_mods_h(data1a))

Fixed Effects

Parameter	AIC	BIC	Chi2	df	p	Converged
m7	-1951	-1749	88.22	56
m7h3	-1904	-1747	157.51	67	< .001
m7h1	-1873	-1729	193.71	70	< .001

Model Performance

rbind(
  performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "All"),
  performance::performance(update(m7_all, data = data1a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1a"),
  performance::performance(update(m7_all, data = data1b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1b"),
  performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(update(m7_all, data = data7a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7a"),
  performance::performance(update(m7_all, data = data7b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7b"),
  performance::performance(update(m7_all, data = data7c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7c"),
  performance::performance(update(m7_all, data = data8a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8a"),
  performance::performance(update(m7_all, data = data8b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 8b"),
  performance::performance(update(m7_all, data = data9), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 9"),
  performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data10a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10a"),
  performance::performance(update(m7_all, data = data10b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 10b")
) |>
  display()

Chi2	RMSEA	CFI	SRMR	Sample
2080.93	0.03	0.98	0.02	All
88.22	0.04	0.98	0.03	Sample 1a
66.50	0.02	0.99	0.03	Sample 1b
123.77	0.04	0.97	0.03	Sample 2
114.03	0.04	0.98	0.03	Sample 3
81.35	0.02	0.99	0.02	Sample 4
83.52	0.05	0.92	0.06	Sample 5
152.88	0.05	0.97	0.03	Sample 6
106.90	0.04	0.97	0.04	Sample 7a
188.21	0.03	0.98	0.02	Sample 7b
122.56	0.04	0.97	0.03	Sample 7c
115.89	0.03	0.98	0.02	Sample 8a
83.61	0.03	0.98	0.03	Sample 8b
2003.31	0.04	0.98	0.02	Sample 9
82.59	0.05	0.97	0.03	Sample 10a
91.47	0.05	0.97	0.04	Sample 10b

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

# Add empirical variables
add_empirical <- function(data, df) {
  x <- data |>
    mutate(
      HungryThirsty = (df$Hungry + df$Thirsty) / 2,
      UrinateDefecate = (df$Urinate + df$Defecate) / 2,
      MusclesPain = (df$Muscles + df$Pain) / 2,
      HeartBreathing = (df$Heart + df$Breathing) / 2,
      ItchBruise = (df$Itch + df$Bruise) / 2,
      WindBurp = (df$Wind + df$Burp) / 2)
  if("Cough" %in% names(df)) {
    mutate(x, CoughSneeze = (df$Cough + df$Sneeze) / 2)
  } else {
    mutate(x, CoughSneeze = df$Sneeze)
  }
}

# Deal with missing data
temp <- predict(m7_all, newdata=data7c)
scores7c <- data.frame(matrix(ncol = 7, nrow = nrow(ega_scores_7c))) 
scores7c[!is.na(ega_scores_7c$EGA1), ] <- temp
names(scores7c) <- names(as.data.frame(temp))

scores <- list(
  sample1a = cbind(ega_scores_1a, data_addprefix(as.data.frame(predict(m7_all, newdata=data1a)), "CFA_")) |>
    add_empirical(data1a) |> 
    mutate(Sample = "Sample 1a"),
  sample1b = cbind(ega_scores_1b, data_addprefix(as.data.frame(predict(m7_all, newdata=data1b)), "CFA_")) |>
    add_empirical(data1b) |> 
    mutate(Sample = "Sample 1b"),
  sample2 = cbind(ega_scores_2, data_addprefix(as.data.frame(predict(m7_all, newdata=data2)), "CFA_")) |>
    add_empirical(data2) |> 
    mutate(Sample = "Sample 2"),
  sample3 = cbind(ega_scores_3, data_addprefix(as.data.frame(predict(m7_all, newdata=data3)), "CFA_")) |>
    add_empirical(data3) |> 
    mutate(Sample = "Sample 3"),
  sample4 = cbind(ega_scores_4, data_addprefix(as.data.frame(predict(m7_all, newdata=data4)), "CFA_")) |>
    add_empirical(data4) |> 
    mutate(Sample = "Sample 4"),
  sample5 = cbind(ega_scores_5, data_addprefix(as.data.frame(predict(m7_all, newdata=data5)), "CFA_")) |>
    add_empirical(data5) |> 
    mutate(Sample = "Sample 5"),
  sample6 = cbind(ega_scores_6, data_addprefix(as.data.frame(predict(m7_all, newdata=data6)), "CFA_")) |>
    add_empirical(data6) |> 
    mutate(Sample = "Sample 6"),
  sample7a = cbind(ega_scores_7a, data_addprefix(as.data.frame(predict(m7_all, newdata=data7a)), "CFA_")) |>
    add_empirical(data7a) |> 
    mutate(Sample = "Sample 7a"),
  sample7b = cbind(ega_scores_7b, data_addprefix(as.data.frame(predict(m7_all, newdata=data7b)), "CFA_")) |>
    add_empirical(data7b) |> 
    mutate(Sample = "Sample 7b"),
  sample7c = cbind(ega_scores_7c, data_addprefix(scores7c, "CFA_")) |> 
    add_empirical(data7c) |> 
    mutate(Sample = "Sample 7c"),
  sample8a = cbind(ega_scores_8a, data_addprefix(as.data.frame(predict(m7_all, newdata=data8a)), "CFA_")) |>
    add_empirical(data8a) |> 
    mutate(Sample = "Sample 8a"),
  sample8b = cbind(ega_scores_8b, data_addprefix(as.data.frame(predict(m7_all, newdata=data8b)), "CFA_")) |>
    add_empirical(data8b) |> 
    mutate(Sample = "Sample 8b"),
  sample9 = cbind(ega_scores_9, data_addprefix(as.data.frame(predict(m7_all, newdata=data9)), "CFA_")) |>
    add_empirical(data9) |> 
    mutate(Sample = "Sample 9"),
  sample10a = cbind(ega_scores_10a, data_addprefix(as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data_all), newdata=data10a)), "CFA_")) |>
    add_empirical(data10a) |> 
    mutate(Sample = "Sample 10a", CFA_CoughSneeze = NA),
  sample10b = cbind(ega_scores_10b, data_addprefix(as.data.frame(predict(m7_all, newdata=data10b)), "CFA_")) |>
    add_empirical(data10b) |> 
    mutate(Sample = "Sample 10b")
)

save(scores, file = "../data/scores.RData")